
TR = function(x.delta, x.delta.d, x.op.y, x.op.d, x.pi, z, d, y, reg.sol, 
              ipw.sol, weights) {
    
    p0.y.est  = reg.sol$p0.y.est
    p0.d.est  = reg.sol$p0.d.est
    f.z.x.est = ipw.sol$f.z.x.est
    
    # DR estimation of beta
    dr.beta.objective = function(beta) {
        obj = t(x.delta.d * weights) %*% ((d - z * tanh(x.delta.d %*% beta) - 
                                     p0.d.est) * 
                                    (2 * z - 1)/f.z.x.est)
        return(sum(obj^2))
    }
    startpars = rep(0, ncol(x.delta.d))
    opt = Optimize(dr.beta.objective, startpars)
    beta.est = opt$par
    delta.d.est = tanh(x.delta.d %*% beta.est)
    
    # DR estimation of alpha
    dr.alpha.objective = function(alpha) {
        obj = t(x.delta * weights) %*% ((y - d * tanh(x.delta %*% alpha) - 
                                             p0.y.est + 
                                   tanh(x.delta %*% alpha) * p0.d.est) * 
                                  (2 * z - 1)/f.z.x.est)
        return(sum(obj^2))
    }
    startpars = rep(0, ncol(x.delta))
    opt = Optimize(dr.alpha.objective, startpars)
    alpha.est = opt$par
    delta.est = tanh(x.delta %*% alpha.est)
    
    # Bounded DR estimation of alpha
    dr.b.alpha.objective = function(alpha) {
        gx = cbind(x.delta, 1/delta.d.est)[,tmp.index]
        obj = t(gx * weights) %*% ((y - d * tanh(x.delta %*% alpha) - p0.y.est + 
                              tanh(x.delta %*% alpha) * p0.d.est) * 
                             (2 * z - 1)/f.z.x.est)
        return(sum(obj^2))
    }
    dim  = ncol(x.delta)+1
    mins = matrix(NA,dim,dim)
    for(i in c(1,dim)){
        startpars = rep(0, ncol(x.delta))
        tmp.index = -i
        opt = Optimize(dr.b.alpha.objective, startpars)
        mins[i,] = c(opt$par,opt$value)
    }
    mins[dim,dim] = mins[dim,dim] + 0.001
    b.alpha.est = mins[which.min(mins[,dim]),1:ncol(x.delta)]
    b.delta.est = tanh(x.delta %*% b.alpha.est)
    
    # TR estimation of Delta
    Delta.est = sum(weights * (  ( y - d*delta.est - p0.y.est + delta.est * p0.d.est) * (2 * z - 1)/ (f.z.x.est * delta.d.est) + delta.est  )) / sum(weights)
    b.Delta.est = sum(weights * b.delta.est) / sum(weights)
    
    return(list(beta.est=beta.est, delta.d.est=delta.d.est,
                alpha.est=alpha.est, delta.est=delta.est, 
                b.alpha.est=b.alpha.est, b.delta.est=b.delta.est,
                Delta.est=Delta.est, b.Delta.est=b.Delta.est,
                opt=opt))
    
}








